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Motivated by recent experiments and a theoretical analysis of the gap equation for the propaga- 
tor of Dirac quasiparticles, we assume that the physics underlying the recently observed removal of 
sublattice and spin degeneracies in graphene in a strong magnetic field is connected with the gen- 
eration of both Dirac masses and spin gaps. The consequences of such a scenario for the existence 
of the gapless edge states with zigzag and armchair boundary conditions are discussed. In the case 
of graphene on a half-plane with a zigzag edge, there are gapless edge states in the spectrum only 
when the spin gap dominates over the mass gap. In the case of an armchair edge, however, the 
existence of the gapless edge states depends on the specific type of mass gaps. 

PACS numbers: 73.43.Cd, 71.70.Di, 81.05.Uw 

I. INTRODUCTION 

A graphite monolayer, or graphene, has become a new exciting topic in physics of two-dimensional electronic 
systems.— i^ 3 - A qualitatively new feature of graphene is that its low-energy quasiparticles are described by a relativistic 
2 + 1-dimensional Dirac theoryi 4 ^^ The spinor structure of the corresponding wave functions is a consequence of the 
honeycomb lattice structure of graphene with two carbon atoms per unit cell. When a magnetic field is applied, 
noninteracting Dirac quasiparticles occupy the Landau levels (LLs) with the following energies: 



E n = ±^2nhv 2 F \eB\/c « ±424yftiy/ B[T] K, (1) 

with n = 0, 1, 2, . . .. Here B is the value of the magnetic field orthogonal to the graphene's plane and vp as 10 6 m/s is 
the Fermi velocity. 

Several anomalous properties of graphene are attributed to the presence of the lowest Landau level (LLL), i.e., the 
n = state in spectrum ((T]), whose energy is independent of the field strength. For example, the anomaly manifests 
itself as the phase shift of ir in the quantum magnetic oscillations of the diagonal conductivity. This phase shift can 
be theoretically understood by using either the semiclassical quantization condition for quasiparticles with a linear 
dispersion,— or a microscopic calculation for both masslcss and massive Dirac fermions.— ^ In the Hall conductivity, 
the anomaly results in an unconventional integer quantum Hall (QH) effect with the plateaus at the filling factors 
v = ±4(n + l/2) i 10 i 1:L i 12 i 13 These and other distinct properties of graphene allow one to unambiguously identify the 
Dirac nature of quasiparticles in experiments . 14 ! 15 

While many unusual properties of graphene can be explained by using the framework of a noninteracting Dirac 
theory, the quasiparticle interactions are not negligible. In fact, they are responsible for the appearance of additional 
QH plateaus with the filling factors v = 0, ±1, ±4 that were first reported in Ref. [TH in the case of sufficiently strong 
magnetic fields, B > 20 T (see also Refs. [lUlUll). 

Recently, we proposed a dynamical mechanism^ which is based on the phenomenon of the magnetic catalysis^ that 
could explain the v = and v = ±1 plateaus in the Hall conductivity of grapheneJ^ The subsequent experiment a 17 i 18 
have revealed several additional features of the v = and v = ±1 plateaus that seem to require modifications of the 
scenario in Ref. l2ll . Among them, the most important is a rather peculiar dissipative nature of the diagonal transport 
at the v — plateau. This seems to suggest that the origin of the v = plateau is associated with a spin gap rather 
than a mass gapi 17 i 23 This conclusion is supported by the fact that the activation energy at the v — plateau is 
vanishingi 17 i 18 Additionally, the diagonal transport is suggested to be dominated by gapless edge states, which should 
exist when the lowest Landau level is split by a large spin gap. 17 ' 23 

Concerning the v = ±1 plateaus, the measurements of the thermal activation energy AE(i> = ±1) point to its 
connection with orbital dynamics. Indeed, the activation energy depends only on the perpendicular component of the 
magnetic ficld j 16 i 18 The dynamical nature of the v = ±1 plateaus is also suggested by the fact that AE(u = ±1) is 
proportional to -\/B . 16 i 18 

Note that, in contrast, the v = ±4 plateaus can be consistently associated with the Zeeman splitting of the n — 1 
Landau level. The corresponding activation energy AE(v = ±4) linearly depends on the total magnetic field and has 
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FIG. 1: Illustration of the lowest Landau level splitting needed to explain v = and v = 1 plateaus in QHE in graphene. 

the same magnitude as the Zeeman energy , 16 ' 18 

Ez = y^bB -0.67 B[T}K, (2) 

where fis — eh/(2mc) is the Bohr magneton and — 2 is the Lande factor in graphene. 

Theoretically, the v = and v = ±1 plateaus come from lifting the approximate degeneracy of the four sublcvels 
at LLL. The degeneracy is a consequence of the "flavor" £7(4) symmetry of the low-energy continuum description of 
graphene in the absence of a Zeeman interaction. This symmetry operates in the space of the sublattice-valley and 
spin degrees of freedom. If it is accepted that the v = plateau is due to a spin gap, then the v = ±1 plateaus should 
result from breaking the sublattice-valley symmetry. This seems to be in agreement with the observations in Ref. HH 

There are essentially two approaches that consider various possibilities of breaking the approximate £7(4) symmetry 
of graphene (for a brief review, see Ref. l24h . 

(i) The quantum Hall ferromagnetism (QHF) scenario 25 ' 26 ' 27 which is connected with the theory of exchange-driven 
spin-splitting of Landau levels in Ref. l28l It exploits an analogy between the four-fold degeneracy of LLs in 
graphene, which is associated with the Z7(4) symmetry, and the 5t/(4) ferromagnetism previously studied in 
the bilayer quantum Hall systems.— In this scenario the QH plateaus with all integer values of the filling factor 
v occur in sufficiently clean samples. The QHF order parameters are described by the densities of conserved 
charges connected with the diagonal generators of the SU(4) C £7(4) symmetry group. 

(ii) The magnetic catalysis (MC) scenario 21 ' 30 ' 31 ' 32 that uses the idea of a spontaneous symmetry breaking due to 
the exciton (chiral) condensation . 22 ' 33 ' 34 ' 35 Such a condensation produces a nonzero Dirac mass term in the 
low-energy theory of graphene. (Note that originally the magnetic catalysis scenario in graphene was motivated 
by the early experiments in highly oriented pyrolytic graphite^) 

As emphasized in Ref. l2lL the plateau v = could appear due to either an enhanced spin gap or a mass term. An 
enhanced spin gap breaks the approximate £7(4) symmetry down to the £7(2)_ x £7(2)+ subgroup which operates in 
the sublattice-valley space and does not mix spin-up (s = +) and spin-down (s = — ) states. A nonzero Dirac mass 
term breaks the symmetry down to another £7 (2) x £7(2)' subgroup, which operates in the spin space. Either of them 
is sufficient to partially lift the four-fold degeneracy of the LLL that is needed in the v = QH state. The structure 
of the energy sublevels at LLL in the case of a nonzero spin gap is illustrated in the left panel of Fig. [T] 

In order to explain the v = ±1 QH plateaus, two different order parameters are required. (Note that the choice 
of two order parameters with given symmetry properties is not unique^ 7 -) This should already be evident from the 
symmetry arguments alone. For example, the simplest possible structure of the energy sublevels for the v — +1 state 
is shown in the right panel of Fig. [TJ The corresponding splitting is possible only if the £7(2)_ x £7(2) + symmetry is 
further reduced, e.g., at least down to the £7(2)_ x £7(1)+ x £7(1)+ subgroup. However, this would not be possible 
without having an additional order parameter that breaks the sublattice-valley symmetry, which is described by the 
simple Lie group SU(2) + C £7(2)+. 

An approach that combines both QHF and MC mechanisms in a unifying scheme was recently proposed in Ref. l37l . 
By making use of a multi-parameter variational ansatz for the quasiparticle propagator, it was found that QHF (fx s 
and p,s) and MC (A s and A s ) order parameters necessarily coexist. In terms of symmetry, the order parameters of 
the first type, i.e., fi s and A s with noncqual values for s = ±, break the £7(4) symmetry down to £7(2)_ x £7(2)+ just 



3 



like the Zeeman term. The order parameters of the other type, fi s and A s , are triplets with respect to the SU(2) S 
group, which is the largest non-Abclian subgroup of the U(2) s . Thus, when either p s or A s has a nonzero vacuum 
expectation value, the symmetry SU(2) S is further broken down to U(l) s . 

The motivation for the present work is to address the question of compatibility of the microscopic dynamics described 
in Ref. [3?] with the gapless edge states, which are apparently needed in the v = state i 17 i 18 Our main results are 
as follows. In the case of graphene on a half-plane with a zigzag edge, there are gapless edge states in the spectrum 
only when the spin gap dominates over the mass gap. In the case of an armchair edge, however, the existence of the 
gapless edge states depends on the specific types of mass gaps. As will be discussed below, these results could have 
important consequences for understanding dynamics in the QH effect in graphene. 

The paper is organized as follows. In Sec. [II] we present a model Lagrangian that captures the most general 
dynamical situation with QHF and MC order parameters, as proposed in Ref. 1371 - The spectrum of the corresponding 
Dirac equation in an external magnetic field is analyzed in Sec. Mil The edge states for zigzag and armchair edges 
are considered in Sees. IIVI andfVl respectively. The main results of the paper are discussed in Sec. IVI1 



II. MODEL WITH DYNAMICAL GAPS 



The low-energy quasiparticle excitations in graphene are described in terms of a four-component Dirac spinor 
^ T s = VipK + As:'>pK + Bs,'4'K-Bs,'4>K-As) ■ The spinor (with a given spin index s = ±) combines the Bloch states on the 
two different sublattices (A and B) of the hexagonal graphene lattice and with the momenta near the two inequivalent 
Dirac points (K+ and K-) of the two-dimensional Brillouin zone. The quadratic part of low-energy Lagrangian density 
for quasiparticles can be written in a relativistic form as 

C = h^ s {t, r) (ij°d t + iv F ^D x + iv Fl 2 D v ) *,(t, r) + £ mass + ^ (p s p s + fi a p t ) , (3) 

s=± s=± 

where = ^7° i s the Dirac conjugated spinor and the operators p s = ^ s j 0l $> s and p s = ^> s j j 5l, Ss are densities 
of conserved charges connected with the chemical potentials p s and p, s (s = ±), respectively. Notice that here the 
Fermi velocity vf ~ c/300 plays the role of the speed of light. The orbital effect of a perpendicular magnetic field 
B = V x A is included via the covariant derivative Di = di + (ie/hc)Ai, where i = x,y and — e < is the electron 
charge. Here, we assume that the vector potential is taken in the Landau gauge: A x = —By and A y = 0, where B is 
the magnitude of the magnetic field. The mass term £ maS s is defined below. 

The 4x4 matrices 7" furnish a reducible representation of the Dirac algebra. Here, we use the following represen- 
tation: 

7° = n <g> t , 7* = -if* ® n, (4) 

where the Pauli matrices f, and Ti (as well as the 2x2 unit matrices fo and To) act on the valley (K+, K_) and the 
sublattice [A, B) indices, respectively. This representation is derived from a tight-binding model for graphene^ It is 
particularly convenient for our purposes in this study because it provides a simple form of the boundary conditions at 
zigzag and armchair edges. As follows from definition (j4]), the 7-matrices satisfy the usual anticommutation relations 
{7^, 7"} = 2<? M1/ , where g^ v = diag(l, —1, —1, —1). Since the matrix 7 5 = i7°7 1 7 2 7 3 is diagonal, 

this representation is conventionally called chiral. Note that the chirality here is identified with the valley index (K + 
or K-)M. 

The general expression for the mass term £ maS s in the Lagrangian density may include singlet (A s ) as well as 
triplet (A s ) contributions with respect to the valley symmetry group SU(2) S . The appearance of the mass term can 
be attributed, for example, to the MC mechanism. In the representation used here, its explicit form reads^ 

Anass = *«(*' r ) ( A ^V - A s7 3 ) *,(*,!■). (6) 

s=± 

Under the time reversal symmetry, the operators associated with the mass parameters A s and A s are odd and even, 
respectively. Concerning the triplet mass term A S 'I' S 7 3 \E' S , it can also be written in other equivalent forms, e.g., as 
A s * s i7 5 ^ s or A s § s * s .ia The latter, in particular, is the usual Dirac mass term. All of these representations are 
equivalent because they are related by the transformations of the SU(2) S symmetry group. For our purposes, however, 
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it is most convenient to use the form in Eq. ([6]) which, as we shall see below, has a simple interpretation in the tight 
binding model. 

In Lagrangian density [Eq. (j^J)], the chemical potentials ji a and jl s allow us to accommodate the QHF order 
parameters in the dynamical model of Ref. l37l . Regarding the transformation properties of /x s and jl s under the flavor 
symmetry, they are similar to those of A s and A s , respectively. 

Before proceeding with further analysis, it is instructive to address the physics interpretation of the mass parameters 
and chemical potentials in more detail. To this end, let us write down the explicit expressions for the corresponding 
operators in the Lagrangian density in terms of separate Bloch components of the spinors as follows: 

A s : *s7 3 *s = 'Wk+aAk+As + tpK_ A8 ^K-As ~ Wk+Bs^K+Bs ~ W K _ Bs 1pK-Bs, (?) 

A s : * s7 3 7 5 * s = ^ K+As ^K+As ~ ll>\c_ As ll>K-Aa ~ ^K+Bs^K+Ba + ^>K_ Bs i>K-Bs, (8) 

jj s : 'J' s 7 7 5 * s = 1p J< K+As ll>K + As ~ 4>k_ As *PK-As + ^K +Bs tpK+Bs ~ ^K_ Bs 1pK-Bs, (9) 

fl s : # s7 °* S = i> f K+As 4>K +As + lpK_ As lpK_Aa + ^K +Bs ^K + Bs + 1pK_ Bs 1pK-Ba- (10) 

Here the operators on the right hand side are linear combinations of the electron densities at specified valleys (K + 
or KJ) and sublattices (A or B). These operators enter into the Lagrangian density together with the parameters 
A s , A s , Us, and p, s , which play the role of Lagrange multipliers. Therefore, the values of the masses and chemical 
potentials control the relative concentrations of electrons at different valleys and sublattices. They are determined 
from the gap equations for Dirac quasiparticles^ 

As seen from Eq. ([7|), the triplet Dirac mass A s is related to the density imbalance between the A and B sublattices. 
Its spontaneous generation leads to a state with a charge density wave i 21 i 30 ' 31 i 32 ' 33 If the values of the dynamical masses 
are nonequal for different spins s = ±, an admixture of an antiferromagnetic wave develops in the ground stated 
Recent scanning tunneling spectroscopy revealed a mass gap near the Dirac point in a single layer graphene sample 
suspended above a graphite substrate4i The gap could be interpreted as a Dirac mass gap induced by a substrate 
perturbation that breaks the sublattice symmetry— In the case of epitaxial graphene grown on SiC, the presence of 
a nonzero Dirac gap is strongly supported by the angle resolved photoemission spectroscopy measurements^ 3 - 

The value of the singlet Dirac mass A s [see Eq. ([8])] controls a mixed density imbalance at the two valleys and the 
two sublattices. Similarly, the chemical potential fl s is connected with the density imbalance between the two valleys 
and, at last, fi s is the usual chemical potential related to the total density of electrons with a given spin. 

In general, in a finite geometry sample, the magnitude of the exchange and Hartree interactions, which determine 
the values of the parameters A s , A s , /x s , and fl s , is likely to vary with the distance from the edges and should be 
calculated in a self-consistent way. The present study of the edge states is done by assuming uniform gaps and uniform 
chemical potentials. We believe, however, that such an idealized treatment should be sufficient to capture the main 
qualitative (although not quantitative) features of the edge states. 



III. DIRAC EQUATION IN AN EXTERNAL MAGNETIC FIELD 



In this section, we study the spectrum of the low-energy quasiparticles in the model of graphene with the most 
general set of parameters A s , A s , fi s , and fl s . The corresponding Dirac equation takes the following form: 



ry u ft<9 t + ih,VFl L D x + ihv F ^D y + ^7° + /i7°7 D + A^y^ - A7 3 $(t, r) = 



..3. ,5 



(11) 



For brevity of notation, the spin index is omitted here and below. For the energy eigenvalue solutions ^(t, r) 
e~ %Et ' h ^{r), the equation reduces to 



hv F {-aiiD x - a 2 iD y ) - fx- fn 5 - iAj 1 -/ 2 + Aa 3 *(r) = ^(r), 



where the a-matrices are 



? 

oti = 7V = 



^ 
ov 



(12) 



(13) 



By using the representation for the 7-matrices in Eq. we can rewrite the Dirac equation in the components as 
follows: 

- A(") -hv F (iD x + D y ) 
-hv F (iD x - D y ) + A ( ~) 

-p(-> - A(+) hv F (iD x + D y )\f ^ BK _ 
Hv F (iD x - Dy) + A(+) )\il> 



IpAK^ 
I^BK^ 



AK- 



= E 
= E 



IpA 



BK4 



IpBK^ 
IpAK. 



(14) 

(15) 
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Here, we introduced the shorthand notation: fj,& = n ± p, and A (±) = A ± A. As we can see, the equations for 
different valleys decouple. This is a very useful property that considerably simplifies the analysis. In each of the two 
decoupled sets of equations, we can express the B-components in terms of the ^-components of the spinors, 



BK+ = 



BK- 



Hv F (iD x - Dy) 
~E + ^+) - A(") 
hv F (iD x + D y ) 



E + /x(")+A(+) 

Then, at K + and K_ valleys, the two-component spinors can be written in the following form: 



ll>K+ = M 



1pAK+ 
hvp(iD x —D y ) I 







Ao 



l/lAK_ 



4>AK_ 



Here, the constants A12 are determined by the normalization conditions, 

d 2 ripK ± (r, k, n)ip K± (r, k', n') = S n ,n'S{k 



k') 



(16) 
(17) 

(18) 
(19) 

(20) 



where k, k' and n, n' are the quantum numbers (e.g., the wave vector along the x or y direction and the Landau level 
index) that characterize the eigenstates of Dirac quasiparticlcs in the magnetic field. 

As follows from Eqs. (fT4f and (fl5|) . the ^-components of the spinors satisfy the following second order differential 
equations: 



(-l 2 D 2 x - l 2 D 2 y + 1) i> AK+ = 2X+^ak+, 
{~l 2 D 2 x - l 2 D 2 y - 1) i> AK _ = 2X^ AK _. 



(21) 



(E + fi^) 2 — (A^) 2 /eq, the Landau energy scale 



Here wc introduced the two dimcnsionless parameters X± 

€q = ^2hvp\eB\/c, and the magnetic length I = y/hc/\eB\. 

In the Landau gauge (A x ,A y ) = (—By,0), the differential equations in Eq. (|2"Tj) do not explicitly depend on the 
x-coordinate, and therefore, the wave functions are plane waves in the x-dircction, 



ipAK + (r, k) = 
iPak„ (r, k) = 



1 



1 



e lhx u + (y,k), 



■>Pbk + 

IpBK- 



1 



Akx 



I2ld 
1 

/2ttT 



v+(y,k), 



e* kx v„(y,k), 



where the functions u±(y,k) depend only on a single combination of the variables. £ 
following equations: 



y/i 



(df-£ 2 Tl + 2A ± ) u± (£) = 0. 
In accordance with Eq. (|18p . the eliminated components v±(y, k) = v±(£) are given by 



v±(0 



(22) 

fcZ, and satisfy the 
(23) 

(24) 



In an infinite system without boundaries, normalizable solutions to Eq. (|23p arc expressed in terms of the Hcrmite 
polynomials, l>(£) oc e - ^ / 2 i7„(£), provided the parameters A± take nonnegative integer values, i.e., 



A± = n, where n = 0, 1, 2 . 



(25) 



Note that the value of the energy E = —/j,( + > + corresponds to a normalizable LLL state in the K + valley. 

For such a state, the apparent singularity in the u+(0 component of the wave function [see Eq. (|24p] is removed 
by a proper redefinition of the normalization constant. The same is not true, however, for the value of the energy 

= state in the K + valley has 
state in the K- valley has 



E 



energy E 
energy E 



( — ) — A^ + ) in the K_ valley. In fact, a direct analysis shows that the only 
^ = — + A' - ) and resides solely on the B sublattice, while the only n 
1^' + A' + ' and resides solely on the A sublattice. 



FIG. 2: Graphene lattice with zigzag and armchair edges. 



IV. EDGE STATES FOR THE ZIGZAG EDGE 



There exist many studies of edge states in graphene under various conditions! 13 ' 17 i 23 ' 39 i 44 ' 45 i 46 ' 47 ' 48 ' 49 ' 50 Here we 
consider a graphene monolayer on the half-plane y > with a zigzag edge parallel to x, as shown in Fig. O To obtain 
the energy spectrum we need to supplement the differential equations for the u±(y,k) and v±{y,k) functions with 
suitable boundary conditions. Such conditions can be derived from the tight-binding modeh 23 i 46 ' 48 For example, for 
a zigzag edge parallel to the x axis, the wave function on the A atoms should vanish at y = 0, 

u+(y = 0) = u_(y = 0) = 0. (26) 

The general solution to Eq. ((23]) is expressed in terms of the parabolic cylinder (Weber) functions U (a, z) and V(a, z),— 

MO-ft B+ ""- At 'V (^ft:,vse) +c ^(i^± V5 S ), (27, 

M«)=Cw(-i±f^) +C4 W-> o+ A<+V (-l±|^). (28, 

Here, for convenience of further analysis, the integration constants Ci and C4 are introduced together with the 
additional factors (E + //+) - A<~)) /e and (i5 + // _ ) + A (+ )) /e , respectively. 

In an infinite system without edges, the normalizable wave functions contain only the parabolic cylinder U(a, z)- 
functions, which are bound at z — > ±00, provided a = — n — 1/2 and n is a nonnegative integer. In this case, 
the following relation is valid: U(—n — 1/2, z) = 2~"/ 2 e~ z / 4 H n (z/\/2), where H n (z) are the Hermite polynomials. 
Therefore, as stated in Sec. IIII1 the spectrum is given by A± = n where n = 0, 1, 2, . . .. (A special nature of LLL 
should be kept in mind: at n = there are only two rather than four possible energy eigenvalues that correspond to 
normalizable states.) 

By using the following recurrent relations for parabolic cylinder functions^ 

U{a,z) = -(a + l)u(a + l,z), 



K dz 2 

and Eq. (|M|) , we obtain the v±(£) functions 



dz 2 J x ' ' 
d z 
dz 2 
d z 
dz' + 2 

^V{a,z) = (a- l -)V{a-l,z), 



U(a,z) = U{a- l,z), 



V(a,z)=V(a + l,z), 



(29) 



mo = -<w (-^ ^) - ^ + ^;; A( "V (-^, ^ ) , (so) 

M0 - C 3 E + ^ e - A(+) U tt^.rt) + <*V (i^, V2, ) . (31) 
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On a half-plane, the normalizable wave functions are also given in terms of only U (a, z)-function, which falls off 
exponentially as z — ■» +00, while the function V(a, z) is growing exponentially in both directions z — > ±00. Therefore, 
we must take C2 = and C\ = 0. In contrast to the case of an infinite plane, on a half-plane, there is no restriction 
for the parameter a to be a negative half-integer. 

With C2 = C4 = 0, the zigzag boundary conditions [ ([26]) ] lead to the following system of equations 



Ci [E + - A<->J D A+ -i(-V^fci) = 0, 
C 3 D X _{-V2kl) = 0. 

Here wc introduced another parabolic cylinder function, D v {z)^ which is related to function U{a, z) in a simple way, 

U(a,z) = D_ a _ y2 (z). (33) 

There are two types of nontrivial solutions that satisfy the boundary conditions f(|32p]. First, by taking C\ 7^ and 
C3 = 0, we find that the equation for the eigenvalues is reduced down to E = — /jW + A' - ) or 

I. Dx+-i(-klV2) = 0. (34) 

The solutions of this type have wave functions with a support only in the K+ valley, 

£ + - AM 

(35) 



= C\ ~ D\ y (v^c) , 

eo v / 

v+(0 = -CxDx. (V2Y) , 



and = «_(£) = 0. The other class of solutions is such that C\ = and C3 7^ 0, and the energy eigenvalues 

satisfy the following equation: 



E + u (_) - A (+) / r- \ 



II. D x _(-kW2) = 0. (36) 
The wave functions for this type of solutions are nonvanishing only in the A'_ valley, i.e., 

II. =C 3 i? A _ (V2c) , 

,(-)_A(+) „ / /r A (37) 

£0 

and = = 0. By making use of the general properties of the parabolic cylinder function D v {z), we can 

understand some qualitative features of the energy spectrum even without solving the equations numerically. To this 
end, we need to know only that, for real v and z, the function D v (z) has no real zeros when v is negative, and has 
exactly [v + 1] real zeros when v is nonnegative^ Here [v + 1] denotes the integer part of v + 1. In view of this 
property, the necessary condition for Eq. (|34[) to be satisfied is A+ > 1. By also including the possibility of the 
dispersionless mode, which is determined by E = —/i^ + A*- - ), wc sec that the complete spectrum in the K + valley 
(solutions of type I) has the following general structure: 

E (k) = V +) +A ( -\ 

I 9 ( 38 ) 

E n (k) = -^ (+) ± y\+(kl,ri)4 + (A(-)) , where \+(kl,n)>l, 

where n = 1, 2, ... is an index that labels different branches of solutions. By making use of the asymptotic behavior 
of the parabolic cylinder functions, one can show that X + (kl, n) ~ n when kl 3> 1. This is expected since large values 
of kl correspond to the states in the bulk, whose wave functions are localized around £ ~ or equivalently y/l ~ kl. 
(In a system without edges, the index n is identified with the usual Landau level index.) 

Similarly, we can constrain the form of the spectrum in the K_ valley (solutions of type II) . The necessary condition 
for having a real solution to Eq. (|36|) is A_ > 0. Thus, the energy spectrum in the K_ valley has the following general 
structure: 



E n {k) = ± \f\-(kl,n)e 2 + (A(+)) 2 , where A_(fcZ,n)>0, (39) 

where n = 0, 1,2, Again, one can show that X-(kl, n) ~ n when kl 3> 1. 
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FIG. 3: The numerical solutions of Eq. (|34|l for the dimensionless parameter A+ (solid line) and Eq. (|36|l for the dimensionless 
parameter A_ (dashed line) in the case of a zigzag boundary. The solid line at A+ = corresponds only to E — — + A' - ' 
solution. 



Our numerical results for X± as functions of kl arc presented in Fig. [3] The solid and dashed lines represent A+ 
and A_, respectively. As expected, there exists an infinite tower of solutions that correspond to an infinite tower of 
Landau levels on a half-plane. In Fig. [3j we show only the first 11 solutions. We also added the constant solution 
A + = that, strictly speaking, represents only the dispersionless mode with the energy E = —fj,( + > + [see the 

first expression in Eq. (j38|) ]. (Formally, A+ = may also mean that E = — — but this is not an energy 

eigenvalue.) 

By analyzing the structure of the spectrum together with the actual dependence of X± on the wave vector, we 
can now determine when gapless modes exist in the spectrum of graphene on a half-plane with a zigzag edge. From 
Eqs. (|55|) and ([55]) , we see that the necessary condition to have a zero energy state is that at least one of the following 
inequalities is satisfied: 

K+ valley: \^ | > yfi* + (A(-)) 2 , (40) 
K_ valley: |//")| > |A (+) |. (41) 

From the fact that there exist branches with A + ~ 1 and A_ ~ at kl 1, we see that this is also the sufficient 
condition. 

An important point to emphasize here is that nonzero masses do not prevent the existence of the gapless edge states 
when the absolute value of A( + ) is less than the absolute value of /x' - ' at least for one choice of the spin. This is 
very similar to the conditions on a graphene ribbon of finite width^ except that there are no edge states associated 
with the second edge in the present work. Our results generalize the findings of previous studies on a half-plane ) 17 i 23 
where only the case with a single nonzero order parameter (either mass or spin gap) was considered. 

Two specific examples of energy spectra, with and without gapless modes, are given in Fig. 3J In the left panel, 
the first few Landau levels in the case of a small spin gap, which is modeled by /i± — =F0.02eo with the subscript 
index denoting the spin, and a larger singlet mass, which is given by A± = ±0.08eo, are shown. Since \^~^ \ < |A' + )|, 
there are no gapless modes in this case. In the right panel of Fig. [4j the low-energy spectrum is shown for another 
choice of parameters, i.e., /i± = =F0.08eo and A± = ±0.02eo, which satisfies the condition in Eq. (|4"Tj) . As expected, 
in this case there are gapless edge states in the spectrum. By taking into account the fact that the group velocities 
of gapless modes, v x = dE/dk\E=o, have opposite signs along the cc-direction, the up- and down-spin states carry 
counter-propagating currents! 17 i 23 It is also curious to note that these gapless states are chiral since they belong to a 
single valley (K-). 

Before concluding this section, it might be appropriate to mention that the examples of spectra shown in Fig.[2]may 
have a direct application to the case of graphene in a strong magnetic field. The corresponding choice of parameters 
with singlet, rather than triplet masses was taken in the same form as in the ground state around the neutral Dirac 
point, which was proposed in the dynamical model of Ref . [37l. In fact, the spectra would look nearly the same also in 
the case of triplet masses, except perhaps for an overall shift of the dispersionless modes, which depend not only on 
the absolute value but also on the sign of the mass terms. 
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FIG. 4: Numerical results for the energy spectra of the first few Landau levels near a zigzag edge of graphene in the case 
of nonzero spin splitting and nonzero singlet masses. The values of parameters are fj,± — =p0.02eo and A± = ±0.08eo in the 
left panel, and fj,± — =p0.08eo and A± = ±0.02eo in the right panel. (The subscript indices in (i± and A± denote the spin 
orientations.) In the first case |/v '| < |A^~'| and there are no gapless modes, in the second case |/u'~~^| > |A' + ^| and gapless 
modes are present. Spin-up and spin-down states are denoted by red (s = +) and blue (s = — ) color of the lines. In the lowest 
energy sublevels the spins are also marked by arrows. The spectra around K+ (K-) point are shown by solid (dashed) lines. 



V. EDGE STATES FOR THE ARMCHAIR EDGE 

In this section, we analyze the spectrum of edge modes in the case of an armchair edge. We take the armchair 
edge parallel to the y-dircction, as shown in Fig. [2] In this case, it is convenient to use a different Landau gauge with 
(A x ,A y ) = (0,Bx). Accordingly, the solutions of Eq. (j2"Tj) are translation invariant along the y-direction, 

iPak + (r, k) = —L=e tky u + (x, k), iPbk + = -^=e lky v+(x, k), 

V2nl V2ttZ (42) 

1>AK-(r,k) = -=e ikv U-(x,k), ip BK _ = — =e lky V- (x, k). 
\Z2ttI \/2irl 

Then, the corresponding differential equations for functions u±(x, k), which are rewritten in terms of the dimensionless 
variable £ = x/l + kl, coincide with Eq. (|23j). The expressions for the eliminated components v±(£), however, slightly 
differ from Eq. (|24|) . and are given by 

" ±(0 = T, ^ + /i w T ^))- (43) 

We consider a graphene sheet in the half-plane x > 0. Since the armchair edge has lattice sites of both A and B 
types, the wave function should vanish at both these sites along the x = line ; 23 ! 46 ' 48 

u + (x = Q)+u_(x = 0) = Q, 

v+{x = 0) +v-(x = 0) = 0. [ ' 

Note that armchair boundary conditions mix the chiralities associated with the K+ and K- valleys. The general 
solutions for the it±(£) functions have the same form as in Eqs. (|27p and (|2"5|) . 

S + M (+)-A(") r /1-2A+ /-A _ /1-2A+ 



«+({) = ft \~ U (^f^ V^J +C 2 V [^f ± , V2£j , (45) 

c / + ^ ) + A(+ V f-I±^, viz) , 



u_(0 - C 3 U ( — , V2£j + C,—^— V I V2£l , (46) 
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FIG. 5: Numerical results for the energy spectra of the first few Landau levels near a armchair edge of graphene in the case 
of nonzero spin splitting and nonzero singlet masses. The values of parameters are fj,± = =p0.02eo and A± = ±0.08eo in the 
left panel, and /U± = =F0.08eo and A± = ±0.02eo in the right panel. (The subscript indices in fj,± and A± denote the spin 
orientations.) In both cases, there are gapless modes in the spectrum. Spin-up and spin-down states are denoted by red (s = +) 
and blue (s = — ) color of the lines. In the lowest energy sublevels the spins are also marked by arrows. 



but with £ = x/l + kl. By using the relations in Eqs. (|43|) and (|29|) . we also obtain the explicit expression for u±(£) 
functions, 

. , Cl U (-1±|^±, ^e) + iC2 E + »<l + *-' v (-i±^±, V5«) , ,47, 

..(0 . iC3 E+ ^-* M v (1^, + icv (1^, Vk) . (48, 

As in the zigzag case, here, normalizable wave functions are given in terms of only the U (a, z)-function, which falls 
off exponentially as z — > +00, unlike the function V(a,z), which grows exponentially in both directions z — > ±00. 
Therefore, we set C2 = and C4 = 0. Then, the armchair boundary conditions [Eq. (|4"4"]) ] lead to the following system 
of equations: 



+ ) r- r- 

C x £ D x+ -i(V2kl) + C 3 Dx + (V2kl) = 0, 

_ R 4- u (-) - A (+) ^ 
C 1 D x _(V2kl) + C 3 — — Da -i(\/2fc/) = 0, 



(49) 



where again we used relation (|33[) to rewrite the expression in terms of the parabolic cylinder function D u (z). This 
system has nontrivial solutions when the determinant of coefficient functions is zero, i.e., 

- A(")) (.E + u(-> - A(+)) / /- \ / \ { r \ f r \ 

± ^ -r>A+-i (V2fcZJ D A _-i (\/2fc^J - L» A+ (V2fc/J D x _ (V2kl) = 0. (50) 

The numerical solutions to this equation for several representative choices of parameters are shown in Figs. O and [5] 

The two cases with singlet masses arc illustrated in Fig. In the left panel, the first few Landau levels in the case of 
\x± = =F0.02eo and A± = ±0.08eo are shown. In the right panel, instead, the corresponding values are /i± = =F0.08eo 
and A± = ±0.02eo. Note that here jl± = A± = 0. (Here, we restored the subscript indices which denote the 
quasiparticle spin orientations.) As we can see, in both cases the spectra contain gapless edge states. This is in strong 
contrast to the zigzag edge case. Indeed, for the armchair edge, gapless modes exist irrespective of the actual relation 
between the values of the singlet masses and spin splitting gaps. In part, this property could be understood from 
the topology of the spectra around the edge and the fact that the singlet mass does not break the SU(2) S valley 
symmetry. The double degenerate sublevels with a given spin, which should exist in the bulk because of the SU(2) S 
symmetry, repel in opposite directions near the edge. Then, gapless modes become almost inevitable at the edge. 

We note that the gapless edge states in Fig. [5] consist of a pair of opposite spin states, carrying counter-propagating 
currents along the edge. This is qualitatively the same situation as found in Ref. [23|. Interestingly, though, if the 
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FIG. 6: Same as in Fig. but for the case of nonzero triplet masses. The values of parameters are (j,± = =p0.02eo and 
A± = 0.08eo in the left panel, and fj,± = =p0.08eo and A± = 0.02eo in the right panel. The existence of gapless modes depends 
on the relative magnitude of \fJ,± \ and |A±|. 
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FIG. 7: Numerical solutions of Eq. (|51[) for the dimensionless parameter A in the case of an armchair boundary. This is valid 
for a general choice of A and fj,, but only if p and A vanish. 



values of singlet masses A + and A_ had the same signs, the opposite spin states would carry currents in the same 
direction along the edge. The observational implications of this fact could be quite unusual. It is not clear, however, 
if such a state can be realized since the dynamical model of Ref. [37] indicates that singlet masses A + and A_ should 
have opposite signs in the ground state. 

The two cases with triplet masses are illustrated in Fig. [6l The values of the parameters in these cases are (i) 
M± = T0.02e and A± = 0.08e (left panel in Fig. EJ) and (ii) //± = =F0.08e and A± = 0.02e (right panel in Fig. EJ. 
These energy spectra resemble the spectra for the zigzag edge, studied in Sec. [TV] There are no gapless edge states 
when the mass is larger than the spin splitting, and there are such states when the mass is smaller than the spin 
splitting. 

In fact, in the case of the triplet mass A and a nonzero \i (but vanishing p and A), we can study the energy spectra 
around the armchair edge in a general case, just like we did for the zigzag edge. In this particular case, the spectral 
equation (|50"|) takes the following simple form: 

\Dl_! (V2kl) - D\ (V2kl\ = 0, (51) 

where A = \(E + /i) 2 — A 2 ]/eg. By expressing A in terms of squares of parabolic cylinder functions from Eq. (|5ip . we 
see that solutions to this equation exist only with A > 0. Therefore, the energy spectrum takes the following form: 



E n {k) = -/i±A/A(/cZ,n)e 2 + A 2 , where A(fcZ,n)>0, (52) 
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where n = 0, 1, 2, Additionally, one can show that X(kl, n) ~ n when |fc^| ^> 1 and k is negative. Our numerical 

results for A as a function of kl are presented in Fig. [7] By combining the numerical information with the general 
expression for the energy (|52|) . we see that the necessary and sufficient condition for having gapless modes is > |A|. 

VI. DISCUSSION 

In this paper, we studied the spectra of edge states in graphene on a half-plane with zigzag and armchair boundary 
conditions, and derived the conditions for the existence of the gapless edge states for various types of masses and 
chemical potentials that could be spontaneously generated in QHE, e.g., at v = and v = ±1 plateaus. 

Our analysis of singlet and triplet Dirac masses [with respect to the valley symmetry group SU(2) S ] shows that 
spectral properties of zigzag and armchair edges arc affected by (i) the relative magnitude of the masses and chemical 
potentials, and (ii) the types of masses. In particular, we found the criteria for the existence of gapless edge states in 
the spectra. These can be summarized as follows. 

(i) Zigzag edge: the necessary and sufficient condition to have a gapless state is that at least one of the following 
inequalities is satisfied: 




(ii) Armchair edge: 

(a) gapless edge states exist always when there are singlet Dirac masses, irrespective of the actual relation 
between the values of the masses and the chemical potentials; 

(b) in the case of triplet Dirac masses, gapless edge states exist if \/i±\ > |A±|, and do not exist otherwise. 

These conditions are consistent with the two limiting cases, analyzed in Ref . |23|. Also, the results in this paper extend 
our previous findings in the case of a graphene ribbon with zigzag edges^ The situation on a half-plane with a zigzag 
edge is essentially the same one, modulo the fact that there is one edge instead of two. 

The above criteria are derived for ideal, smooth edges and for a perfect graphene layer without disorder. In reality, 
the available graphene samples are disordered. Because of the geometrical roughness and impurities, they do not 
have perfect zigzag or armchair edges either. Then, the corresponding boundary conditions for the graphene wave 
functions may be different from those used in the current study^ Additionally, the bonds of the carbon atoms at the 
edges can be saturated by foreign atoms modifying even perfectly smooth and regular edges^ Therefore, it is of great 
importance to study the effects of various types of disorder in graphene. This is, however, beyond the scope of the 
present paper. Here we limit our study to an idealized model in order to provide a clean benchmark calculation before 
a more detailed investigation of disorder is undertaken. By taking into account a considerable improvement in sample 
quality seen in graphene suspended above a graphite substrate^ or above a Si/Si02 gate electrode^, it is possible 
that the clean limit already provides a reasonable qualitative description of edge states. Additionally, because of the 
special nature of the LLL, the role of some types of disorder may be strongly suppressed^ For example, the effect 
of the randomness in the bond couplings and in the on-site potential caused by short range interactions is studied in 
Ref. [56|. It is shown that the degeneracy of K± points is not lifted by the on-site disorder, but can be removed by the 
randomness in the bond couplings. 

The results here are of interest in connection with the interpretation of the v = Hall plateau. Indeed, the gapless 
edge states should play an important role in the charge transport of graphene in a strong magnetic field. Their presence 
is expected to make graphene a so-called quantum Hall metal, while their absence should make it an insulator i 17 i 23 
The actual temperature dependence of the longitudinal resistivity at the v = plateau in Refs . Il6| and Il7l is consistent 
with the metal type. This conclusion may be disputed in view of the recent data from Ref. |2(J that reveal a clear 
plateau at v = 0, but the temperature dependence of the diagonal component of the resistivity signals a crossover to 
an insulating state in high fields. The latter observations do not seem to support the existence of gapless edge states. 

Our analysis in this paper as well as in Ref. [33 suggests that the conditions for the existence and absence of gapless 
edge states sensitively depend on the values of QHF and MC order parameters that characterize the nature of the 
corresponding QH state. Moreover, the microscopic analysis of Ref. [37] indicates that the order parameters of both 
types necessarily coexist. Therefore, the dynamics is very likely to be rich and full of surprises. The situation with 
the edge states is probably just one of such surprises. 



13 



Acknowledgments 

The authors acknowledge useful discussions with E.V. Gorbar, H. Fertig, I.F. Herbut, M.I. Katsnelson, L. Levitov, 
and B.I. Shklovskii. V.P.G. and S.G.S. thank A.K. Gcim for the discussion of the experimental data that indicate the 
existence of the gapless edge states in graphene. The work of V.P.G. was supported by the SCOPES Project No. IB 
7320-110848 of the NSF-CH, Grant No. 10/07-N "Nanostructure systems, nanomaterials, nanotechnologics," and the 
Program of Fundamental Research of the Physics and Astronomy Division of the National Academy of Sciences of 
Ukraine. The work of V.A.M. was supported by the Natural Sciences and Engineering Research Council of Canada. 



* On leave from Bogolyubov Institute for Theoretical Physics, 03680, Kiev, Ukraine 

1 K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, and A.A. Firsov, Science 
306, 666 (2004). 

2 A.K. Geim and K.S. Novoselov, Nat. Mater. 6, 183 (2007). 

3 A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov and A.K. Geim, larXiv:0709.TT63l Rev. Mod. Phys. (to be 
published) . 

4 G.W. Semenoff, Phys. Rev. Lett. 53, 2449 (1984). 

5 D.P. DiVincenzo and E.J. Mele, Phys. Rev. B 29, 1685 (1984). 

6 F.D.M. Haldane, Phys. Rev. Lett. 61, 2015 (1988). 

7 G.P. Mikitik and Yu.V. Sharlai, Phys. Rev. Lett. 82, 2147 (1999). 

8 S.G. Sharapov, V.P. Gusynin, and H. Beck, Phys. Rev. B 69, 075104 (2004). 

9 V.P. Gusynin and S.G. Sharapov, Phys. Rev. B 71, 125124 (2005). 

10 Y. Zheng and T. Ando, Phys. Rev. B 65, 245420 (2002). 

11 V.P. Gusynin and S.G. Sharapov, Phys. Rev. Lett. 95, 146801 (2005). 

12 V.P. Gusynin and S.G. Sharapov, Phys. Rev. B 73, 245411 (2006). 

13 N.M.R. Peres, F. Guinea, and A.H. Castro Neto, Phys. Rev. B 73, 125411 (2006). 

14 K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, and A.A. Firsov, 
Nature (London) 438, 197 (2005). 

15 Y. Zhang, Y.-W. Tan, H.L. Stormer, and P. Kim, Nature (London) 438, 201 (2005). 

16 Y. Zhang, Z. Jiang, J.P. Small, M.S. Purewal, Y.-W. Tan, M. Fazlollahi, JD. Chudow, J. A. Jaszczak, H.L. Stormer, and P. 
Kim, Phys. Rev. Lett. 96, 136806 (2006). 

17 DA. Abanin, K.S. Novoselov, U. Zeitler, PA. Lee, A.K. Geim, and L.S. Levitov, Phys. Rev. Lett. 98, 196806 (2007). 

18 Z. Jiang, Y. Zhang, H.L. Stormer, and P. Kim, Phys. Rev. Lett. 99, 106802 (2007). 

19 A.J.M. Giesbers, U. Zeitler, M.I. Katsnelson, LA. Ponomarenko, T.M. Mohiuddin, and J.C. Maan, Phys. Rev. Lett. 99, 
206803 (2007). 

20 J.C. Checkelsky, L. Li, and N.P. Ong, Phys. Rev. Lett. 100, 206801 (2008). 

21 V.P. Gusynin, VA. Miransky, S.G. Sharapov, and LA. Shovkovy, Phys. Rev. B 74, 195429 (2006). 

22 V.P. Gusynin, VA. Miransky, and LA. Shovkovy, Phys. Rev. Lett. 73, 3499 (1994); Phys. Rev. D 52, 4718 (1995); Nucl. 
Phys. B 462, 249 (1996). 

23 DA. Abanin, PA. Lee, and L.S. Levitov, Phys. Rev. Lett. 96, 176803 (2006); Solid State Commun. 143, 77 (2007). 

24 K. Yang, Solid State Commun. 143, 27 (2007). 

25 K. Nomura and A. H. MacDonald, Phys. Rev. Lett. 96, 256602 (2006). 

26 M.O. Goerbig, R. Moessner, and B. Doucot, Phys. Rev. B 74, 161407(R) (2006). 

27 J. Alicea and M.P.A. Fisher, Phys. Rev. B 74, 075422 (2006). 

28 M.M. Fogler and B.I. Shklovskii, Phys. Rev. B 52, 17366 (1995). 

29 In both cases, the S'[/(4) symmetry is associated with spin-pseudospin degrees of freedom, but in graphene, pseudospin is 
related to the sublattice index, in bilayer quantum Hall systems, pseudospin is related to the layer index. 

30 I.F. Herbut, Phys. Rev. Lett. 97, 146401 (2006); I.F. Herbut, Phys. Rev. B 75, 165411 (2007); I.F. Herbut, ibid. 76, 085432 
(2007). 

31 J.-N. Fuchs and P. Lederer, Phys. Rev. Lett. 98, 016803 (2007). 

32 M. Ezawa, J. Phys.Soc. Jpn. 76, 094701 (2007); Physica E (Amsterdam) 40, 269 (2007). 

33 D.V. Khveshchenko, Phys. Rev. Lett. 87, 206401 (2001); D.V. Khveshchenko, ibid. 87, 246802 (2001); D.V. Khveshchenko 
and H. Leal, Nucl. Phys. B 687, 323 (2004); D.V. Khveshchenko and WF. Shively, Phys. Rev. B 73, 115104 (2006). 

34 E.V. Gorbar, V.P. Gusynin, VA. Miransky, and LA. Shovkovy, Phys. Rev. B 66, 045108 (2002). 

35 E.V. Gorbar, V.P. Gusynin, VA. Miransky, and LA. Shovkovy, Phys. Lett. A 313, 472 (2003). 

36 Y. Kopelevich, V.V. Lemanov, S. Moehlecke, and J.H.S. Torres, Fiz. Tverd. Tela 41, 2135 (1999) [Phys. Solid State 41, 1959 
(1999)]; H. Kempa, Y. Kopelevich, F. Mrowka, A. Setzer, J. H. S. Torres, R. Hohne, and P. Esquinazi, Solid State Commun. 
115, 539 (2000); M.S. Sercheli, Y. Kopelevich, R. R. da Silva, J.H.S. Torres, and C. Rettori, Solid State Commun. 121, 579 
(2002); Y. Kopelevich, J. C. Medina Pantoja, R. R . da Silva, F. Mro wka, and P. Esquinazi, Phys. Lett. A 355, 233 (2006). 

37 E. V. Gorbar, V. P. Gusynin, and V. A. Miransky. larXiv:0710.3527l (unpublished) . 



14 



38 V.P. Gusynin, S.G. Sharapov, and J. P. Carbotte, Int. J. Mod. P hys. B 21, 4611 (2 007). 

39 V.P. Gusynin, V.A. Miransky, S.G. Sharapov, and LA. Shovkovv. larXiv:0801. 07081 (unpublished - ) . 

40 Strictly speaking, in order to preserve the SU(2) S valley symmetry of the model, all three components of the triplet should 
be included in the Lagrangian density on equal footing. The ground state, however, will correspond to a specific "vacuum 
alignment", e.g., characterized by a nonvanishing vacuum expectation value of the operator *l! a^ 3 s . 

41 G. Li, A. Luican, and E.Y. Andrei. larXiv:0803.4016l (unpublished') . 

42 G. Giovannetti, P. A. Khomyakov, G. Brocks, P.J. Kelly, and J. van den Brink, Phys. Rev B 76, 073103 (2007). 

43 S.Y. Zhou, G.-H. Gweon, A.V. Fedorov, P.N. First, W.A. de Heer, D.-H. Lee, F. Guinea, A.H. Castro Neto, and A. Lanzara, 
Nat. Mater. 6, 770 (2007). 

44 K. Nakada, M. Fujita, G. Dresselhaus and M.S. Dresselhaus, Phys. Rev. B 54, 17954 (1996). 

45 M. Fujita, K. Wakabayashi, K. Nakada and K. Kusakabe, J. Phys. Soc. Japan 65, 1920 (1996). 

46 E. McCann and V.I. Fal'ko, J.Phys.: Condens. Matter 16, 2371 (2004). 

47 C.L. Kane and E.J. Mele, Phys. Rev. Lett. 95, 146802 (2005); ibid. 95, 226801 (2005). 

48 L. Brey and H. A. Fertig, Phys. Rev. B73, 235411 (2006); ibid. 73, 195408 (2006). 

49 N.M.R. Peres, A.H. Castro Neto and F. Guinea, Phys. Rev. B 73, 241403(R) (2006). 

50 J. M. Pereira, F.M. Peeters and P. Vasilopoulos, Phys. Rev. B 75, 125433 (2007). 

51 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, 
(U. S. GPO, Washington, D.C., 1972), p. 685. 

52 H. Bateman and A. Erdelyi, Higher Transcendental Functions (Mc Graw-Hill, New York, 1953), Vol. II. 

53 I. Martin and Ya.M. Blanter, arXiv:0705.0532 A.R. Akhmerov and C.W.J. Beenakker, Phys. Rev. B 77, 085423 (2008); 
V. Cvetkovic and Z. Te sanovic, arXiv:08 0TT2T2l 

54 S. Dutta and S.K. Pati. larXiv:0712.4130l (unpublished). 

55 K.I. Bolotin, K.J. Sikes, Z. Jiang, G. Fudenberg, J. Hone, P. Kim, and H. L. Stormer. lar"Xiv:0802.2389l (unpublished). 

56 M. Koshino and T. Ando, Phys. Rev. B 75, 033412 (2007). 



